4 


! 

NASA Technical Memorandum 101473 
ICOMP-89-3 


Numerical Computation of Shock Wave- 
Turbulent Boundary Layer Interaction in 
Transonic Flow Over an Axisymmetric 
Curved Hill 


S.-W. Kim 

Institute for Computational Mechanics in Propulsion 
Lewis Research Center 
Cleveland, Ohio 

(MSA-TM-1C1473) 80UEE1CAL CCt GfclAlIOB Cf tj89-i11S2 

SECCK KAVE-IOBEBIEKT ECOKIAti IllEB 
1 F1EBAC1 ICh lb IGitiSCKlC fIXS CUi IS 

AlISYEMElitlC CUFVIE Hlli. ( FA £ A ) 3 1 p Unclas 

CSCi 200 G3/34 0191958 


February 1989 






NUMERICAL COMPUTATION OF SHOCK WAVE - TURBULENT BOUNDARY LAYER INTERACTION 
IN TRANSONIC FLOW OVER AN AXISYMMETRIC CURVED HILL 


S.-W. Kim* 

Institute for Computational Mechanics in Propulsion 
Lewis Research Center 
Cleveland , Ohio 44135 

Summary 

A control -volume based finite difference computation of a turbulent 
transonic flow over an axisymmetric curved hill is presented. The numerical 
method is based on the SIMPLE algorithm, and hence the conservation of mass 
equation is replaced by a pressure correction equation for compressible 
flows. The turbulence is described by a k-e turbulence model supplemented by 
a near-wall turbulence model. In the method, the dissipation rate in the 
region very close to the wall is obtained from an algebraic equation and that 
for the rest of the flow domain is obtained by solving a partial differential 
equation for the dissipation rate. The other flow equations are integrated up 
to the wall. It is shown that the present turbulence model yields the correct 
location of the compression shock. The other computational results are also 
in good agreement with experimental data. 


Work funded under Space Act Agreement C99066G. 



Nomenclature 


coefficient for axial velocity correction equation 

coefficient for radial velocity correction equation 

constant coefficient for f^ equation (-0.025) 

constant coefficient for f^ equation (-0.00001) 

constant coefficient for f € equation 

turbulence model constants for e equation, (i-1,2) 

constant coefficient for eddy viscosity equation (“0.09) 

wall damping function for eddy viscosity equation 

wall damping function for c w equation 

turbulent kinetic energy 

effective thermal conductivity (=k m +kt : ) 

thermal conductivity 

turbulent thermal conductivity (=Cp^ t /o*p) 
pressure 

production rate of turbulent kinetic energy 
gas constant. 

turbulent Reynolds number (=k^/(i/e i ) ) 

radial coordinate 

temperature 

Reynolds averaged velocity in axial direction 
friction velocity (=7(r w /p)) 

Reynolds stress 
velocity vector (-(u,v)) 

Reynolds averaged velocity in radial directions 
axial coordinate 


normal distance from the wall 
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wall coordinate (**u r y/i/) 
dissipation rate 

dissipation rate of turbulent kinetic energy 

dissipation rate inside the near-wall layer 

von Karman constant (=0.41) 

molecular viscosity 

effective viscosity (=/i+/i t ) 

turbulent viscosity 

kinematic viscosity of fluid 

turbulent eddy viscosity 

density 

turbulent Prandtl number for k-equation 

turbulent Prandtl number for energy equation 

turbulent Prandtl number for e -equation 

wall shearing stress 

dissipation function for energy equation 


Superscripts 

* current value 

' incremental (or corrective) value 

A non-dimensional value normalized by the free stream value 

Mathematical symbol 

X summation 
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Introduction 


The transonic flow over an axisymmetric curved hill [1] has received 
considerable attention in recent years as a bench mark test case to assess 
the capability of numerical methods as well as turbulence models to be used 
as design/analysis tools for fluid machinery. In the numerical calculation 
of turbulent transonic flows, prediction of the correct pressure field 
depends on the location of the shock and the location of the shock depends 
on the viscous force development. Therefore, emphasis was laid upon 
selecting and/or developing a suitable turbulence model in References 2 and 
3. The turbulence models tested in these references ranged from algebraic 
turbulence models to two-equation turbulence models, and varying degrees of 
success have been obtained. 

The transonic flow calculation presented herein constitutes one of the 
earliest applications of the newly developed numerical method [4] as well 
as the turbulence model [5]. These are described below. 

The control -volume method based on the SIMPLE algorithm [6,7] is 
mostly used to solve incompressible flows the domain of which can be 
discretized by an orthogonal mesh. Due to their strongly convergent nature, 
pressure correction methods have been used extensively to solve complex 
turbulent flows including chemically reacting turbulent flows [8-9]. The 
numerical method used herein is an extension of the pressure correction 
method to solve incompressible as well as compressible flows with 
arbitrary, complex geometries. The compressible flow equations are mostly 
solved by flux splitting methods. The Beam-Warming method [10] and the 
McCormack method [11] are representatives of the flux splitting methods. 

The flux splitting methods have originally been developed to solve the 
Euler equations and then extended to include the viscous term to solve the 
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Navier-Stokes equations. The most distinguished difference between the two 
classes of methods lies in the way the diffusion term is treated. In the 
pressure correction methods, the diffusion term is incorporated into the 
stiffness matrix while, in the flux splitting methods, the diffusion term 
is incorporated into the system of equations as the load vector term. For 
turbulent flows with extensive recirculation zones, the pressure correction 
methods may be numerically more stable than the other class of methods, 
conceptually; however, the pressure correction methods have mostly been 
used for incompressible flows and the flux splitting methods have mostly 
been used for compressible flows. Therefore, definitive advantages and 
disadvantages of these two classes of methods can not be discussed with 
confidence as yet. 

A few papers to extend the SIMPLE method to solve compressible flows 
have appeared in recent years [4,12-14]. Some difficulties have been 
encountered in the course of these studies. A fully staggered grid layout 
was used in the original SIMPLE method to solve flow equations using 
orthogonal meshes [6,7]. In many flow problems of practical importance, the 
boundary geometries are complex and arbitrary shaped blockages may exist 
inside the flow path. One difficulty was identifying a suitable grid layout 
to solve the Navier-Stokes equations defined on complex geometries. In 
Reference 12, a collocated grid layout was used and an artificial 
dissipation was included to prevent velocity-pressure decoupling. In 
Reference 15, a fully staggered grid layout for incompressible flows was 
used and the velocity vector was located at all grid points except at the 
pressure grid point, hence the number of degrees of freedom for velocities 
is doubled and that of pressure remains the same as in the original case. 

In References 4 and 14, the velocities were located at the same grid points 


5 



and the pressure was located at the centroid of the four adjacent velocity 
grid points. This grid layout has been used successfully in the penalty 
finite element method for a long time [16]. It was first used in the 
control -volume based finite difference method by Vanka et. al. [17]. They 
mentioned that it was not ea ;y to obtain convergent solutions due to 
velocity-pressure decoupling. The mechanism that leads to the 
velocity-pressure decoupled solution was heuris tically shown in Reference 
15. In Reference 14, the velocity-pressure decoupling was eliminated by 
using a non- conforming domain for mass imbalance calculation. In Reference 
4, the velocity-pressure decoupling was eliminated by treating the pressure 
correction equation as a standard partial differential equation rather than 
treating it as a constraint condition. In the method, the off-diagonal 
terms were moved to the load vector terra and the resulting system of 
equations was solved using the tri-diagonal matrix algorithm (TDMA) . Thus 
any uncertainty that may arise due to the use of a non-conforming domain 
for mass imbalance calculation does not exist in the present method [4], 
Another difficulty was identifying the most suitable numerical 
procedure for pressure correction. The SIMPLE-R [6], the SIMPLE-C [18], and 
the standard SIMPLE algorithm [6] were used in References 13, 14, and 4, 
respectively. The pressure, velocity, and density were corrected based on 
the incremental pressure (or pressure correction) in References 13 and 14. 
In Reference 4, only the pressure and velocity were corrected from the 
incremental pressure. Density was obtained from the equation of state for a 
perfect gas so that the same numerical procedure could be equally 
applicable for numerical computation of chemically reacting turbulent flows 
in the future. The numerical procedure for pressure correction, especially 
for compressible flows, may need to be further studied in the future. 
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The accuracy and the convergence nature of the numerical method used 
herein has been tested by solving a class of example flow cases. The 
example problems considered in Reference 4 include: a developing channel 
flow, a developing pipe flow, a two-dimensional laminar flow in a 90 degree 
bent channel, polar cavity flows, and a turbulent supersonic flow over a 
compression ramp. It was found that the numerical method used herein 
yielded accurate computational results even when highly skewed, unequally 
spaced, curved grids were used. It is also found that the present method 
was strongly convergent for high Reynolds number flows as well as for flows 
with complex geometries. 

In numerical calculations of turbulent flows, the near-wall region has 
been incorporated into the numerical analyses usually by using the wall 
functions [19], two- or multi-layer turbulence models [20-22], and low 
Reynolds number turbulence models [23], 

The most widely used wall function methods have been derived from the 
logarithmic velocity profile based on the experimental observation that the 
turbulence in the near-wall region can be described in terms of the wall 
shearing stress. The wall function methods are not valid if the logarithmic 
velocity profile no longer prevails in the near-wall region. Such cases 
include separated and/or reattaching flows, unsteady flows, flows over 
surfaces with mass injection and/or suction, and near-wall low turbulent 
Reynolds number flows. 

In the two- or multi-layer turbulence models, the turbulence in the 
near-wall layer has been described by algebraic equations. In some of these 
turbulence models [20,21], the turbulent kinetic energy and the dissipation 
rate in the near-wall layer have been constructed by piecewise continuous 
functions. As a consequence, the computational results may depend on the 
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location of the partition between the near-wall layer and the fully 
turbulent region. In Reference 20, a quadratic variation of the turbulent 
kinetic energy has been assumed in the near-wall region and the dissipation 
rate has been derived from the logarithmic velocity profile equation. The 
underlying assumption that the turbulent kinetic energy grows in proportion 
to the square of the distance from the wall in the near-wall region is also 
questionable, since the quadratic variation of the turbulent kinetic energy 
is valid only inside the viscous sublayer and becomes invalid as the fully 
turbulent region is approached. 

In the low Reynolds number turbulence models, the high Reynolds number 
turbulence models have been extended to include the near-wall low 
turbulence effect [23]. In this class of methods, the near-wall low 
turbulence effects have been incorporated into the high Reynolds number 
turbulence models by using wall damping functions. These wall damping 
functions have been derived mostly from numerical experiments in such a way 
that the low Reynolds number turbulence models could approximately 
reproduce the experimentally observed turbulent flow phenomena in the 
near-wall region. In this class of methods, a significant number of grid 
points has to be assigned in the near wall region in order to resolve the 
stiff dissipation rate equation. 

Aside from the above two classes of methods, a new approach has been 
used in Chen and Patel [24] to resolve the near-wall turbulence. In the 
method, only the turbulent kinetic energy equation of the k-e turbulence 
model has been extended to include the near-wall low turbulence region and 
the dissipation rate inside the near-wall layer has been prescribed 
algebraically. The dissipation rate equation has been obtained from a 
k-equation turbulence model [23). Thus the turbulent kinetic energy and the 
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dissipation rate vary smoothly from the wall toward the outside fully 
turbulent region. In this case, it would be more appropriate to classify 
this turbulence model as a partially low Reynolds number turbulence model 
to distinguish it from the other two classes of methods. Advantages of the 
partially low Reynolds number approach over the other methods can be 
summarized as follows. The turbulence in external flows and that in 
near-wall boundary layer flows have quite different length scales [26]. The 
turbulence length scale of the external flows is related to the flow field 
characteristics. On the other hand, the turbulence length scale of boundary 
layer flows is strongly related to the normal distance from the wall. This 
characteristic of the wall bounded turbulent flows can be described quite 
naturally by the partially low Reynolds number turbulence models. The same 
purpose could be achieved by the low Reynolds number turbulence models as 
more experimental data become available; however, the gradient of the 
dissipation rate is the most stiff in the near-wall region, and a great 
number of grid points has to be used in this region for the low Reynolds 
number turbulence models to resolve the dissipation rate. Therefore, the 
partially low Reynolds number turbulence models seem to be more 
advantageous as compared with the low Reynolds number turbulence models, 
unless the low Reynolds number turbulence models can describe the wall 
bounded turbulent flows more accurately. 

The turbulence model used in this report belongs to the partially low 
Reynolds number turbulence models. Development of the near-wall turbulence 
model and its application to fully developed turbulent channel and pipe 
flows can be found in Reference 5. It has been shown in the reference that 
the present near-wall turbulence model can resolve the over-shoot phenomena 
of the turbulent kinetic energy and the dissipation rate in the region very 
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close to the wall and that significantly improved computational results for 
the turbulence structure in the near-wall region have been obtained. More 
application of the near-wall turbulence model for complex turbulent flows 
such as the supersonic flow over a compression ramp with shock wave - 
turbulent boundary layer interaction [27] can be found in Reference 4. 

Reynolds Averaged Navier- Stokes Equations and Numerical Method 
The compressible turbulent flow equations are given as; 

a ia 

— (pu) + (prv) - 0. (1) 

dx r dr 
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— (puu) + (prvu) 
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dx r dr 

dx r dr 

dx 
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and the density is obtained from the perfect gas law given as p*=pRT. A 
turbulent Prandtl number (oj) of 0.75 was used for the energy equation, see 
the Nomenclature. The molecular viscosity and the thermal conductivity were 
obtained from the Sutherland's laws given as [28]; 


T ' 

3/2 

T 0 + S' 

To. 


T + S, 


(5) 


where /i Q - 1.716 x 10' 5 Kg/m-sec, T Q * 273.1° Kelvin, S - 110.6° Kelvin; 
and 



f— > 

T 

3/2 

T 0 + S' 

k o 

T 

L x oJ 


T + S, 
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where k Q - 0.0264 Kg/m-K, T Q = 273.1° Kelvin, and S = 194.4° Kelvin. 


In the present method, all flow variables, except pressure, have been 
located at the same grid points and the pressure node has been located at 
the centroid of the cell. The control volume for the pressure correction 
equation is defined as the cell enclosed by the four neighboring grid 
points. Note that in the control -volume based finite difference methods, 
the discrete system of equations is derived by integrating the governing 
differential equations over the control volume [6]. For curvilinear grids, 
the required number of interpolations to obtain flow variables at the cell 
boundaries is significantly reduced by using the present grid layout. 
Enhanced convergence rate is partly attributed to the grid layout which 
required fewer interpolations [4] . 

The pressure correction equation for compressible flows is briefly 


discussed below for completeness. As 
method, the density, the velocities, 

in the standard pressure correction 
and the pressure are decomposed as; 

ll 

* 

+ 

“C> 

(7) 

* . * 

u - u + u' , 

(8) 

V = V* + v' , 

(9) 

p - p* + p' 

(10) 


The incremental pressure is related to the incremental density and the 
incremental velocities as; 

P r * P'RT (11) 
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( 12 ) 


dp' 



dp' 

v' - - Av (13) 

dy 


where eq. (11) has been obtained from the equation of state; and eqs. (12) 
and (13) have been obtained from the discrete u- and v-momentum equations, 
respectively [4,6], Substituting eqs. (7-13) into (1) yields, after some 


rearrangement ; 
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(14) 


where the last term in eq. (14) represents the mass imbalance. The 
incremental pressure is obtained by solving eq. (14), and the corresponding 
incremental velocities are obtained from eqs. (12-13). The flow variables 
are updated using eqs. (8-10), and these updated flow variables are used in 
computing the new current flow variables by solving eqs. (2-4) together 
with the turbulence equations. The iterative solution process is repeated 
until the mass imbalance in eq. (14) becomes negligible. In solving eq. 

(14), the off-diagonal terms in the discrete pressure correction equation 
were moved to the load vector term, and the resulting five diagonal system 
of equations was solved using the TDMA [4,6]. Even the slightest symptom of 
the velocity-pressure decoupling was not observed in the present flow case 
as well as for all the flow cases considered in Reference 4. 
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Tu rbul once Equations 


A k-c turbulence model supplemented with a near-wall turbulence model 
is described below. The turbulent kinetic energy equation for the entire 
flow domain is given as; 


8 Id 

a 

' ak' 

i a 

ak' 

— (puk) + (prvk) 

- — 

^e- 

+ 

r ^e — 

3x r 3r 

dx 

ax 

r dr 

arj 


where the production rate of turbulent kinetic energy (P r ) is the same as 
the dissipation function for the energy equation ($) . 

The dissipation rate inside the near-wall layer is given as; 


£ 1 



where 


W 3 /4 k 3/2 


*y 


f e - 1- exp(-A £ R t ) 


Rf- = 


C..-V2 


A, = 




2^ 


(17) 


Note that in eq.(17) represents the standard dissipation rate for 
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near-wall turbulent flows in equilibrium state. The dissipation rate given 
as eq. (16) is used for eq . (15) in the near-wall region. For y~0 , eq. (16) 

takes the limit value given as 2j/k/y 2 . Slightly away from the wall where 
the turbulence is in the equilibrium state, f e takes the unit value and eq . 
(16) becomes identical to eq. (17). The dissipation rate given as eq. (16) 
is formally identical to the one proposed by Wolfshtein [25]. The 
qualitatively similar dissipation rates used in the k-equation turbulence 
models of Gibson et. al . [29], Hassid and Poreh [30], and Acharya and 

Reynolds [31] can be written as; 


e w “ ^eh 


k 3 / 2 2i/k 


(18) 


where f e ^ is a wall damping function varying from zero on the wall to unity 
in the fully turbulent region. Note that the second term on the right hand 
side of eq. (18) can be obtained as an analytical solution of the turbulent 
kinetic energy equation for a limiting case as y approaches the wall. 
However, this term persists far out into the fully turbulent region, and 
hence eq . (16) has been preferred over eq (18), see Reference 5 for more 

discussion . 

The dissipation rate for the rest of the flow domain is obtained by 
solving the convec tion- dif fus ion equation for the dissipation rate equation 
given as; 

a l a a ( ae] l a ( a«] ? r e e 2 

— (put) + (prve ) = — p e — + rp e — + pci - pc 2 — (19) 

3x r dr ax[ 3xJ r 3r[ dr) k k 
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The turbulence model constants used are given as: <7^=0. 75, o c ~1.15, 
c^=1.39, and C2“1.88. These turbulence model constants approximately 
satisfy the near-wall equilibrium turbulence condition and the decay rate 
of the grid turbulence observed in the experiment [32]. Further discussion 
on the establishment of these turbulence model constants can be found in 
References 33-34. 

The eddy viscosity equation inside the near-wall layer is given as; 
k 2 

V Z " c /if f M (20) 

f l 


where f^-l-exp(A^7R t + A 2 R t ^) . The wall damping function is a linear 
function of the distance from the wall in the viscous sublayer region and 
become unity in the fully turbulent region. The eddy viscosity given as eq. 
(20) grows in proportion to the cubic power of the distance from the wall. 
It can be found in References 5 and 23 that the near-wall analysis yields 
the same growth rate of the eddy viscosity in the region very close to the 
wall. The eddy viscosity in the rest of the flow domain is given as; 



c 


The domain for each differential equation is shown schematically in 
Figure 1 for clarity. For wall bounded turbulent flows, the equilibrium 
region extends from y + «30 up to y + ~300 . Thus the partition between the 
near-wall region and the fully turbulent outer region can be located 
between greater than 100 and less than 300 approximately. 
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Computational Results 

The transonic flow considered in the present study is schematically 
shown in Figure 2. In the experiment [2], an axisymmetric circular-arc bump 
of thickness 1.9 cm and a chord length of 20.3 cm was attached at 60 cm 
downstream of a circular cylinder with the external diameter of 15.2 cm. 

The free stream Mach number was 0.875, the Reynolds number was 13.2xl0^/m, 
and the boundary layer thickness of the approaching transonic flow was 0.01 
meters . 

In the following calculations, the inlet boundary was located at one 
chord length upstream of the forward corner of the bump; and the exit 
boundary, at one chord length downstream of the rear end of the bump. Some 
degree of uncertainty that may be caused by numerical diffusion and 
inadequate grid size can always exist in any numerical analysis. To reduce 
the numerical uncertainty, a relatively fine grid (78 X 53) and a highly 
fine grid (108 X 65) have been used in the present study. The computational 
results obtained using these two grids differed by no more than a few 
percent. The computational results presented herein were obtained using the 
highly fine grid shown in Figure 3. The inlet boundary condition for the 
axial velocity and the turbulent kinetic energy were obtained from 
experimental data for a fully developed flat plate flow [35]. The 
non-dimensional velocity and the turbulent kinetic energy profiles were 
scaled to yield a boundary layer thickness of 0.01 meters at the inlet 
boundary. Uniform static pressure and uniform enthalpy were also prescribed 
at the inlet boundary. The no-slip boundary condition for velocities, a 
vanishing turbulent kinetic energy, and a constant temperature which 
corresponds to the free stream stagnation temperature were prescribed at 


17 


the solid wall boundary. The free stream flow condition was prescribed at 
the top boundary, and a vanishing gradient boundary condition was used for 
all flow variables at the exit boundary. The partition between the 
near-wall layer and the external region was located at approximately 5 
percent of the boundary layer thickness away from the wall and 10 grid 
points were allocated inside the near-wall layer. The initial guess was 
obtained by extending the inlet boundary condition in the axial direction. 

The discrete finite difference system of equations was solved 
iteratively until the error norms became smaller than the prescribed 
convergence criteria. Each iteration consisted of 7 sweeps of the pressure 
correction equation and 3 sweeps for the rest of the flow equations in the 
axial and in the radial directions, respectively. The pressure was updated 
using an under -relaxation factor of 0.57; and the rest of the flow 
variables, using an under- relaxation factor of 0.47. With the use of these 
under-relaxation parameters, divergence or convergence to an erroneous 
solution has not been encountered. The convergence criteria used were; 

N A A 

Rl = £ C | V -(PV)| C < ei (22) 

c=l 

R 2 = |(a? +1 - a. n )/A n+1 | < e 2 , j-l,N, (23) 

i.j i.j i 

where N c is the number of pressure control volumes; the superscript n 
denotes the iteration level; the subscript i={u, v, p, T, k, e} denotes 
each flow variable; the subscript j denotes each grid point; N denotes the 
number of degrees of freedom for each flow variable; and Aj_ denotes the 
maximum magnitude of the i-th flow variable. The iteration was terminated 
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when either eq. (22) or eq. (23) was satisfied. The converged solution was 
obtained after 750 iterations for e^=5xl0"^ and 62 = 1 x 10 *^. At the time of 
convergence and R 2 were approximately 4.5x10"^ and l.lxl0"\ 
respectively. The required computational time was approximately 11 minutes 
for the CRAY-XMP at NASA/LeRC . 

The calculated iso-Mach lines are shown in Figure 4, where the 
incremental Mach number between the contour lines is 0.05. The static 
pressure contour lines are shown in Figure 5, where the pressure has been 
normalized by the inlet total pressure and the incremental pressure between 
the contour lines is 0.01. It can be seen from Figures 4 and 5 that the 
flow field is characterized by a relatively small supersonic pocket 
attached at the top of the curved hill. Neither the iso-Mach lines nor the 
pressure contour are available in References 1-3, and direct comparison of 
these contour lines with other computational results could not be made. The 
calculated static pressure on the wall is compared with experimental data 
as well as the numerical results of References 2 and 3 in Figure 6. It can 
be seen in the figure that the pressure distribution on the wall obtained 
by the McCormack method using Johnson's algebraic turbulence model [3] 
(hereafter abbreviated as J-A model) compares most favorably with 
experimental data. It is shown in the figure that the present turbulence 
model yielded significantly improved pressure distribution over the one [2] 
obtained by the McCormack method using the Wilcox-Rubesin Turbulence model 
(hereafter abbreviated as W-R model) . It can also be seen in the Figure 
that the present turbulence model and the W-R model yielded a slightly more 
spread-out shock than the experimental data. In the present study, the 
computational results are compared mostly with those obtained using the W-R 
model [2] for the following reasons. Firstly, the k-e class turbulence 
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models have wider applicability than the algebraic turbulence models, 
especially for turbulent flows through complex geometries. Secondly, the 
k>£ turbulence models seldom contain a flow* dependent adjustable constant; 
however, this is not always true for algebraic turbulence models [3], And 
lastly, more extensive computational results for the present flow case, 
such as the flow separation and the turbulent kinetic energy profiles, were 
presented in Reference [2] than in Reference [3]. 

The streamline contour at the rear end of the hill is shown in Figure 
7. The measure flow separation zone extended from x/c=0.7 to x/c=l . 1 , where 
c is the chord length. The present method yielded the flow recirculation 
zone extending from x/c=0.81 to x/c=1.12. The recirculation zone obtained 
in Reference 2 using the Cebeci-Smith turbulence model and the W-R model 
extended from x/c=0.81 to x/c-1 . 1 and from x/c=0.81 to x/c=l . 2 , 
respectively. It can be seen that these computational results exhibited 
decent comparison with the experimental data. 

The mean velocity profiles at four axial locations are compared with 
experimental data as well as with those of Reference 2 in Figure 8. It can 
be seen that both computational results exhibited fair comparison with the 
experimental data. At x/c=0.75, the present velocity profile compared more 
favorably with experimental data than that of Reference 2. The velocity 
profiles obtained using the J-A model [3] also compared favorably with 
experimental data. The level of agreement with experimental data was the 
same as the present one. The improvement in the calculated velocity profile 
is attributed to the turbulence models which yielded better turbulence 
fields than the W-R model. 

The calculated turbulent kinetic energy profiles are compared with 
experimental data in Figure 9. It can be seen in the figure that the 
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magnitude and location of the turbulent kinetic energy compared favorably 
with experimental data. Again, the present turbulence model yielded 
significantly improved computational result at x/c=0.75. It has been shown 
in Reference 5 that the near-wall turbulence model can resolve details of 
the near-wall turbulence including the over-shoot phenomena of the 
turbulent kinetic energy and the dissipation rate in the region very close 
to the wall. Considering this fact, it is not fortuitous that the near-wall 
turbulence model yielded improved turbulent kinetic energy profile at 
x/c-0.75. 

The Reynolds stress profiles at the same axial locations are shown in 
Figure 10. It can be seen that the present turbulence model under predicted 
the magnitude of the Reynolds stress. The Reynolds stress profiles 
presented in Reference 3 also exhibited the same trend as the present 
results. It has long been argued that the algebraic turbulence models and 
the k-c turbulence models can not resolve recirculating turbulent flows 
accurately. The same argument can be applied to the present computational 
results. However, the flow field in the downstream region of the curved 
hill did not affect the upstream region significantly and the correct 
surface pressure distribution and the shock location were obtained in the 
numerical calculations . 


Conclusions 

A control -volume based finite difference computation of a turbulent 
transonic flow using a k-e turbulence model supplemented with a near-wall 
turbulence model has been presented. It has been shown that the present 
method was strongly convergent for the high Reynolds number flow with 
arbitrary geometry and that the method yields significantly accurate 
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computational results for the transonic flow with the shock wave -turbulent 
boundary layer interaction. The strongly convergent nature is attributed to 
the near-wall turbulence model, the pressure - staggered grid layout, and the 
pressure correction method. The improved computational results are 
attributed to the near-wall turbulence model. 

In the near-wall turbulence model, the turbulent kinetic energy 
equation was integrated up to the wall and the dissipation rate inside the 
near-wall region was obtained from an algebraic equation. This approach was 
found to be more advantageous than the low Reynolds number turbulence 
models since the stiff dissipation rate equation in the near-wall region 
need not be solved numerically. In some numerical calculations of turbulent 
flows, algebraic turbulence models have been preferred over k-€ turbulence 
models to overcome the numerical difficulty originating from the stiff 
dissipation rate equation [3,36,37], This difficulty is considerably 
reduced by using the partially low Reynolds number turbulence models. More 
importantly, the near-wall turbulence model yielded significantly accurate 
computational results for the near-wall turbulence field of the shock wave 
- turbulent boundary layer interaction flow. For the separated flow region 
at the rear end of the hill, the present turbulence model as well as the 
J-A model yielded degenerated computational results. The degenerated 
computational results are attributed to the limited capability of the 
algebraic turbulence model and the k-c turbulence model to resolve 
recirculating turbulent flows [8,35]. However, the usefulness of the 
present turbulence model to solve turbulent flows with a small separated 
flow region has been demonstrated in this study. 
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FIGURE 1. - DOMAIN FOR EACH PARTIAL DIFFERENTIAL EQUATION. 
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FIGURE 2. - TRANSONIC FLOW OVER AN AX I SYMMETRIC CURVED HILL. 
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FIGURE 3. - DISCRETIZATION OF THE FLOW DOMAIN (108 X 65 MESH). 



FIGURE A. - ISO-MACH LINES. 



FIGURE 5. - PRESSURE CONTOUR. 
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